%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main Exercises
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
clc

% Load parameters
load param

x_s=myx;
[err, TimeSeries, mytargets, mymodel0] = obj_func( x_s ); 

% Calculate moments
gini = cell2mat(TimeSeries(1:4,2));
dissim = cell2mat(TimeSeries(1:4,1));
rankrank = cell2mat(TimeSeries(1:4,6));
n_size = cell2mat(TimeSeries(1:4,7));
ed_returns = cell2mat(TimeSeries(1:4,8));
returns = cell2mat(TimeSeries(1:4,9));
house_price = cell2mat(TimeSeries(1:4,10));
n_income = cell2mat(TimeSeries(1:4,3));
col_share = cell2mat(TimeSeries(1:4,16));
dissm_edu = cell2mat(TimeSeries(1:4,17));
share_rich = cell2mat(TimeSeries(1:4,14));
w_ratio = cell2mat(TimeSeries(1:4,11));
ability = cell2mat(TimeSeries(1:4,20));

tmp1 = [house_price(:,1)./house_price(:,2), house_price(:,2)./house_price(:,3)];
tmp2 = [n_income(:,1)./n_income(:,2), n_income(:,2)./n_income(:,3)];

output = [gini dissim rankrank returns w_ratio share_rich n_size  tmp1  ed_returns col_share dissm_edu ];

%% Figure 8
gini_data = [0.3763; 0.4112; 0.4411; 0.4661];
dissim_data = [0.3341891; 0.3926; 0.4148; 0.4434701];


fig=figure;
hold on
plot([1980 1990 2000 2010], gini,'- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], gini_data,'-- .','LineWidth',1.5, 'MarkerSize',20)
xticks([1980 1990 2000 2010])
%ylim([0.34 0.48])
lgd=legend({'Model','Data'}, 'Location','southoutside','FontSize',12,'Box','off');
lgd.NumColumns = 2;
set(gca,'XGrid','off','YGrid','on')
ax = gca;
ax.FontSize = 12; 
ylabel('Gini')
hold off
set(fig,'Units','Inches');
pos = get(fig,'Position');
set(fig,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(fig,'figures/fig-baselinea.pdf','-dpdf','-r0')

fig=figure;
hold on
plot([1980 1990 2000 2010], dissim,'- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], dissim_data,'-- .','LineWidth',1.5, 'MarkerSize',20)
xticks([1980 1990 2000 2010])
%ylim([0.25 0.45])
lgd=legend({'Model','Data'}, 'Location','southoutside','FontSize',12,'Box','off');
lgd.NumColumns = 2;
set(gca,'XGrid','off','YGrid','on')
ax = gca;
ax.FontSize = 12;
ylabel('Dissimilarity')
hold off
set(fig,'Units','Inches');
pos = get(fig,'Position');
set(fig,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(fig,'figures/fig-baselineb.pdf','-dpdf','-r0')

%% Tables 

Ss = cell2mat(TimeSeries(:,13));
SS_A_B = (Ss(:,1)./Ss(:,2))';  
SS_B_C = (Ss(:,2)./Ss(:,3))'; 

% Table 3
table_3 = [ed_returns'; ...
gini'; ...
dissim'; ...
SS_A_B(1:end-1) ; ... 
SS_B_C(1:end-1) ; ...
tmp1'; ...
n_size']

% Table 5 
table_5 = 100*(cell2mat(TimeSeries(1,21)) + cell2mat(TimeSeries(2,21)) + cell2mat(TimeSeries(3,21)))/3


% Q5-Q5 over time:
Qt_1980 = cell2mat(TimeSeries(1,21));
Qt_1980(5,5)

Qt_2010 = cell2mat(TimeSeries(4,21));
Qt_2010(5,5)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Counterfactuals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x_c = [x_s(1:end-2) 1.75*3-x_s(end-1)-x_s(end) x_s(end-1:end)];

random_mixing    = counter_factual(x_c, 1); % Neighborhood choice probabilites are uniform
fixed_location_1 = counter_factual(x_c, 2); % Fix steady-state neighborhood
fixed_location_2 = counter_factual(x_c, 3); % Fix steady-state neighborhood, prices adjust
fix_loc_edu      = counter_factual(x_c, 4); % Fix location and education to steady state

fixed_spillover     =  decomposition( x_c, 1); % Fix steady-state spillover
mechanical_effect   =  decomposition( x_c, 2); % Fix rents and choices, but spillovers adjust
partial_equilibrium =  decomposition( x_c, 3); % Fix rents, but spillovers and choice adjust
fixed_edu           =  decomposition( x_c, 4); % Fix education to steady state

% No spillover:
x_b = x_s;
x_b(8) = 0;
[err, beta_1_0, mytargets, mymodel0] = obj_func( x_b ); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gini_random = [gini(1); cell2mat(random_mixing(1:3,2))];
gini_fixed  = [gini(1); cell2mat(fixed_location_1(1:3,2))];
gini_fix_loc_edu  = [gini(1); cell2mat(fix_loc_edu(1:3,2)) ];
gini_noedu  = [ cell2mat(fixed_edu(1:4,2)) ];

dissim_random = [dissim(1); cell2mat(random_mixing(1:3,1))];
dissim_fixed  = [dissim(1); cell2mat(fixed_location_1(1:3,1))];
dissim_fix_loc_edu  = [dissim(1); cell2mat(fix_loc_edu(1:3,1)) ];
dissim_noedu  = [ cell2mat(fixed_edu(1:4,1)) ];


gin_no_spill = [ cell2mat(beta_1_0(1:4,2))];
gin_fix_spill  = [cell2mat(fixed_spillover(1:4,2))];

dissim_no_spill = [ cell2mat(beta_1_0(1:4,1))];
dissim_fix_spill  = [cell2mat(fixed_spillover(1:4,1))];

% Figure 9
fig=figure;
hold on
plot([1980 1990 2000 2010], gini,'- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], gini_random,'-- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], gini_fixed,': .','LineWidth',1.5, 'MarkerSize',20,'color',[118,173,72]/255)
xticks([1980 1990 2000 2010])
%ylim([0.34 0.48])
lgd=legend({'Baseline','Random Mixing','Fixed Location'}, 'Location','southoutside','FontSize',12,'Box','off');
lgd.NumColumns = 3;
set(gca,'XGrid','off','YGrid','on')
ax = gca;
ax.FontSize = 12; 
ylabel('Gini')
hold off
set(fig,'Units','Inches');
pos = get(fig,'Position');
set(fig,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(fig,'figures/fig-counter-gini.pdf','-dpdf','-r0')

contr_rand_gini  = 1-(gini_random(end)-gini_random(1))/(gini(end)-gini(1));
contr_fixed_gini = 1-(gini_fixed(end)-gini_fixed(1))/(gini(end)-gini(1));

rankrank_random = [rankrank(1); cell2mat(random_mixing(1:3,6))];
rankrank_fixed = [rankrank(1); cell2mat(fixed_location_1(1:3,6))];
contr_rand_rankrank = 1-(rankrank_random(end)-rankrank_random(1))/(rankrank(end)-rankrank(1));
contr_fixed_rankrank = 1-(rankrank_fixed(end)-rankrank_fixed(1))/(rankrank(end)-rankrank(1));

fig=figure;
hold on
plot([1980 1990 2000 2010], dissim,'- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], dissim_random,'-- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], dissim_fixed,': .','LineWidth',1.5, 'MarkerSize',20,'color',[118,173,72]/255)
xticks([1980 1990 2000 2010])
%ylim([0.34 0.48])
lgd=legend({'Baseline','Random Mixing','Fixed Location'}, 'Location','southoutside','FontSize',12,'Box','off');
lgd.NumColumns = 3;
set(gca,'XGrid','off','YGrid','on')
ax = gca;
ax.FontSize = 12; 
ylabel('Dissimilarity')
hold off
set(fig,'Units','Inches');
pos = get(fig,'Position');
set(fig,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(fig,'figures/fig-counter-segregation.pdf','-dpdf','-r0')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 10
contr_fixspill_gini  = 1-(gin_fix_spill(end)-gin_fix_spill(1))/(gini(end)-gini(1));
contr_fixspill_dissim  = 1-(dissim_fix_spill(end)-dissim_fix_spill(1))/(dissim(end)-dissim(1));

fig=figure;
hold on
plot([1980 1990 2000 2010], gini,'- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], gin_fix_spill,'-- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], gin_no_spill,': .','LineWidth',1.5, 'MarkerSize',20,'color',[118,173,72]/255)
xticks([1980 1990 2000 2010])
%ylim([0.34 0.48])
lgd=legend({'Baseline','Fixed Spillover','No Spillover'}, 'Location','southoutside','FontSize',12,'Box','off');
lgd.NumColumns = 3;
set(gca,'XGrid','off','YGrid','on')
ax = gca;
ax.FontSize = 12; 
ylabel('Gini')
hold off
set(fig,'Units','Inches');
pos = get(fig,'Position');
set(fig,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(fig,'figures/fig-deco-a.pdf','-dpdf','-r0')

fig=figure;
hold on
plot([1980 1990 2000 2010], dissim,'- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], dissim_fix_spill,'-- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], dissim_no_spill,': .','LineWidth',1.5, 'MarkerSize',20,'color',[118,173,72]/255)
xticks([1980 1990 2000 2010])
%ylim([0.34 0.48])
lgd=legend({'Baseline','Fixed Spillover','No Spillover'}, 'Location','southoutside','FontSize',12,'Box','off');
lgd.NumColumns = 3;
set(gca,'XGrid','off','YGrid','on')
ax = gca;
ax.FontSize = 12; 
ylabel('Dissimilarity')
hold off
set(fig,'Units','Inches');
pos = get(fig,'Position');
set(fig,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(fig,'figures/fig-deco-b.pdf','-dpdf','-r0')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 11
Ss_mech = cell2mat(mechanical_effect(:,13));
Ss_pe = cell2mat(partial_equilibrium(:,13));

SS_A_B_mech = Ss_mech(:,1)./Ss_mech(:,2);
SS_B_C_mech = Ss_mech(:,2)./Ss_mech(:,3);

SS_A_B_pe = Ss_pe(:,1)./Ss_pe(:,2);
SS_B_C_pe = Ss_pe(:,2)./Ss_pe(:,3);

contr_ge_spill_ab  = 1-(SS_A_B_pe(4)-SS_A_B_pe(1))/(SS_A_B(4)-SS_A_B(1));
contr_ge_spill_bc  = 1-(SS_B_C_pe(4)-SS_B_C_pe(1))/(SS_B_C(4)-SS_B_C(1));

fig=figure;
hold on
plot([1980 1990 2000 2010], SS_A_B(1:4),'- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], SS_A_B_pe(1:4),'-- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], SS_A_B_mech(1:4),': .','LineWidth',1.5, 'MarkerSize',20,'color',[118,173,72]/255)
xticks([1980 1990 2000 2010])
%ylim([0.34 0.48])
lgd=legend({'Baseline','Partial Equilibrium','Mechanical Effect'}, 'Location','southoutside','FontSize',12,'Box','off');
lgd.NumColumns = 3;
set(gca,'XGrid','off','YGrid','on')
ax = gca;
ax.FontSize = 12; 
ylabel('S_A/S_B')
hold off
set(fig,'Units','Inches');
pos = get(fig,'Position');
set(fig,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(fig,'figures/fig-decomposition-a.pdf','-dpdf','-r0')

fig=figure;
hold on
plot([1980 1990 2000 2010], SS_B_C(1:4),'- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], SS_B_C_pe(1:4),'-- .','LineWidth',1.5, 'MarkerSize',20)
plot([1980 1990 2000 2010], SS_B_C_mech(1:4),': .','LineWidth',1.5, 'MarkerSize',20,'color',[118,173,72]/255)
xticks([1980 1990 2000 2010])
%ylim([0.34 0.48])
lgd=legend({'Baseline','Partial Equilibrium','Mechanical Effect'}, 'Location','southoutside','FontSize',12,'Box','off');
lgd.NumColumns = 3;
set(gca,'XGrid','off','YGrid','on')
ax = gca;
ax.FontSize = 12; 
ylabel('S_B/S_C')
hold off
set(fig,'Units','Inches');
pos = get(fig,'Position');
set(fig,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(fig,'figures/fig-decomposition-b.pdf','-dpdf','-r0')

